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Abstract. We consider numerical approximations of the Monge- Ampere equation 

det D^u = f, f > with Dirichlct boundary conditions on a convex bounded 
domain f2 in M", n = 2, 3. We make a comparative study of three existing methods 
suitable for finite element computations. We construct conforming approximations 
in the framework of the spline clement method where constraints and interelement 
continuities are enforced using Lagrange multipliers. 



This paper addresses the numerical solution of the Dirichlet problem for the Monge- 
Ampere equation 



/ > 0. The domain Q C M"" is a convex domain with polygonal boundary dQ. 

The above equation is a fully nonlinear equation in the sense that it is nonlinear 
in the highest order derivatives. Fully nonlinear equations have in general multi- 
ple solutions, and even if the domain is smooth, the solution may not be smooth. 
For the Monge-Ampere equation, the notion of generalized solution in the sense of 
Alcxandrov-Bakelman and that of viscosity solution [39] arc the best known to give 
a meaning to the second derivatives even when the solution is not smooth. To a 
continuous convex function, one associates the so-called Monge-Ampere measure and 
(1.1) is said to have a solution in the sense of Alexandrov if the density of that mea- 
sure with respect to the Lebesgue measure is equal to /. Continuous convex viscosity 
solutions are defined " in terms of certain inequalities holding wherever the graph 
of the solution is touched on one side or the other by a smooth test function " [44] . 
In the case of (1.1) the two notions are equivalent for / continuous, [39]. We will 
assume throughout this paper that / is continuous on Q and g continuous on dfl. 
Equation (1.1) then has at most two solutions when n = 2, [24] p. 324. and a unique 
generalized solution in the class of convex functions, [1, 23]. In general, the theory of 
viscosity solutions [25, 19, 39] provides a framework for existence and uniqueness of 
solutions of fully nonlinear equations. 

The more general Monge-Ampere equation has form 



1. Introduction 



(1.1) 



det D^u = f in fl, u = g on dfl, 




is the Hessian of u and /, g are given functions with 



(1.2) 



det D^u — H{x, u, Du) in Q, u — g on dQ, 
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where Du denotes the gradient of u and is is a given Hamiltonian, at least con- 
tinuous and nondecreasing in u. They appear in various geometric and variational 
problems, e.g. the Monge-Kantorovich problem, and in kinetic theory. They also 
appear in applied fields such as meteorology, fiuid mechanics, nonlinear elasticity, 
antenna design, material sciences and mathematical finance. A huge amount of liter- 
ature on theoretical questions about these equations is available. A selection in the 
areas cited above include [2, 58, 12, 3, 19, 49, 38, 26, 21]. 

Researchers working on the Monge-Kantorovich Problem, MKP, c.f. [37] for back- 
ground, have noted the problematic lack of good numerical solvers for the Monge- 
Ampere type equations. Following [29], we quote from [15], "It follows from this 
theoretical result that a natural computational solution of the MKP is the nu- 
merical resolution of the Monge- Ampere equation" ..." Unfortunately, this fully non- 
linear second-order elliptic equation has not received much attention from numerical 
analysts and, to the best of our knowledge, there is no efficient finite-difference or 
finite-element methods, comparable to those developed for linear second-order elliptic 
equations (such as fast Poisson solvers, multigrid methods, preconditioned conjugate 
gradient methods,. • • )." 

Existing numerical work on the Monge-Ampere type equations included [50, 42, 20] 
where the generalized solution in the sense of Alexandrov-Bakelman is approximated 
directly. Other works with proven convergence results is [48, 35] where finite difference 
schemes satisfying conditions for convergence of [14] were constructed. There have 
been an explosion of recent numerical results for the Monge-Ampere equations. We 
have the recent papers [45, 46, 59, 18, 52] which do not address adequately the 
situations where the Monge-Ampere equation does not have a smooth solution, c.f. 
Test 2 and Test 3 in Section 4. For progress in this direction we refer to the series of 
papers [28, 29, 27] and the vanishing moment method in [32, 33, 34]. Finite difference 
methods which computes viscosity solutions of the Monge-Ampere equation and an 
iterative method amenable to finite element computations were reported in [16, 36]. 
See also [17] for an optimization approach. 

In this paper, we use the spline element method to compute numerical solution of 
the Monge-Ampere equation. It is a conforming finite element implementation with 
Lagrange multipliers. We will obtain conforming approximations for the three di- 
mensional Monge-Ampere equation. We extend the convergence analysis of Newton's 
method due to [45] to bounded smooth domains using Schauder estimates proved 
in [56]. However [45] did not address the convergence of the discrete approxima- 
tions. We give error estimates for conforming approximations of a smooth solution. 
The Monge-Ampere equation leads to a non-coercive variational problem, a difficulty 
which is partially handled by the vanishing moment method (the parameter e cannot 
be taken equal to 0). We show that for smooth solutions, the spline element method 
is robust for the associated singular perturbation problem. The numerical results 
mainly examine the performance of three numerical methods, Newton's method, the 
vanishing moment method and the Benamou-Proese-Oberman iterative method, on 
three test functions suggested in [27]: a smooth radial solution, a non-smooth solu- 
tion for which no exact formula is known and a solution not in H'^{fl). In this paper. 
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we will refer to the Benamou-Froese-Oberman iterative method as the BFO iterative 
method, which we extend to three dimensions. 

The paper is organized as follows: In the first section, we review the spline ele- 
ment discretization. The following section is devoted to the variational formulations 
associated to Newton's method, for which we give convergence results for smooth 
solutions, and the vanishing moment method. Here we introduce the three dimen- 
sional version of the BFO iterative method. The last section is devoted to numerical 
experiments. We will use C for a generic constant but will index specific constants. 



The sphne element method has been described in [4, 7, 8, 13, 41] under different 
names and more recently in [6]. It can be described as a conforming finite element 
implementation with Lagrange multipliers. We first outline the main steps of the 
method, discuss its advantages and possible disadvantage. We then give more details 
of this approach but refer to the above references for explicit formulas. 

First, start with a representation of a piecewise discontinuous polynomial as a vector 
in M.^ , for some integer > 0. Then express boundary conditions and constraints 
including global continuity or smoothness conditions as linear relations. In our work, 
we use the Bernstein basis representation, [4, 6] which is very convenient to express 
smoothness conditions and very popular in computer aided geometric design. Hence 
the term "spline" in the name of the method. Splines are piecewise polynomials with 
smoothness properties. One then write a discrete version of the equation along with 
a discrete version of the spaces of trial and test functions. The boundary conditions 
and constraints are enforced using Lagrange multipliers. We are lead to saddle point 
problems which are solved by an augmented Lagrangian algorithm (sequences of linear 
equations with size N x N). The approach here should be contrasted with other 
approaches where Lagrange multipliers are introduced before discretization, i.e. in 
[9] or the discontinuous Galerkin methods. 

The spline element method, stands out as a robust, flexible, efficient and accurate 
method. It can be applied to a wide range of PDEs in science and engineering in 
both two and three dimensions; constraints and smoothness are enforced exactly 
and there is no need to implement basis functions with the required properties; it is 
particularly suitable for fourth order PDEs; no inf-sup condition are needed to ap- 
proximate Lagrange multipliers which arise due to the constraints, e.g. the pressure 
term in the Navier-Stokes equations; one gets in a single implementation approxima- 
tions of variable order. Other advantages of the method include the flexibility of using 
polynomials of different degrees on different elements [41], the facility of implement- 
ing boundary conditions and the simplicity of a posteriori error estimates since the 
method is conforming for many problems. A possible disadvantage of this approach 
is the high number of degrees of freedom and the need to solve saddle point problems. 

Let T be a conforming partition of Q into triangles or tetrahedra. We consider a 
general variational problem: Find u eW such that 



2. Spline element discretization 



(2.1) 




for all V eV, 
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where W and V are respectively the space of trial and test functions. We will assume 
that the form / is bounded and linear and a is a continuous mapping in some sense 
onWxV which is linear in the argument v. 

Let Wh and Vh be conforming subspaces of W and V respectively. We can write 

Wh = {ce R^, Rc = G}, Vh = {ce R^, Rc = 0}, 

for a suitable vector G and R a suitable matrix which encodes the constraints on 
the solution, e.g. smoothness and boundary conditions. Here /i is a discretization 
parameter which controls the size of the elements in the partition. 

The condition a{u, v) — {I, v) for all v eV translates to 

K{c)d = L^d yd e Vh, that is for all d with Rd = 0, 

for a suitable matrix K{c) which depends on c and L is a vector of coefficients asso- 
ciated to the linear form I. If for example {l,v) = j^fv, then L^d — d^MF where 
M is a mass matrix and F a vector of coefficients associated to the spline interpolant 
of /. In the linear case K{c) can be written c^K. 

Introducing a Lagrange multiplier A, the functional 

K(c)d-L'^d + X^Rd, 
vanishes identically on Vh- The stronger condition 

K(c) + X^R = L^, 

along with the side condition Rc — G are the discrete equations to be solved. 

By a slight abuse of notation, after linearization by Newton's method, the above 
nonlinear equation leads to solving systems of type 

c^K + X^R = L^. 

The approximation c oi u E W thus is a hmit of a sequence of solutions of systems 
of type 





R^' 




c 




" L ' 


R 







A 




G 



It is therefore enough to consider the hnear case. If we assume for simphcity that V — 
W and that the form a is bilinear, symmetric, continuous and ^-elliptic, existence of a 
discrete solution follows from Lax-Milgram lemma. On the other hand, the ellipticity 
assures uniqueness of the component c which can be retrieved by a least squares 
solution of the above system [4] . The Lagrange multiplier A may not be unique. To 
avoid systems of large size, a variant of the augmented Lagrangian algorithm is used. 
For this, we consider the sequence of problems 

(-) ( -% ) 

where A'-^^ is a suitable initial guess for example A*^°^ = 0, M is a suitable matrix and 
/X > is a small parameter taken in practice in the order of 10~^. It is possible to 
solve for c*-^"'"^^ in terms of c*^'^ A uniform convergence rate in n for this algorithm 
was shown in [5]. 







A(^+i) 





L 



G - /xMA(') 
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3. Variational formulations 

The BFO iterative method for solving (1.1) has a clear variational formulation as it 
consists in solving a sequence of Poisson problems: 

(3.1) Auk+i = y/(Auky + 2(f-detD^Uk), 

for the two dimensional case [16] and extended here to three dimensions 

(3.2) Atx,+i = {{Aukf + 9(/ - dot D'uk))K 
with Uk+i = g on dfl. Since 

(3.3) detD^-u < 4^(Axi)", 

the above formula enforces partial convexity since Ait > is a necessary condition 
for the Hessian of u to be semi positive definite. Wc note that the constant 2 in (3.1) 
may be changed to 4. We next discuss Newton's method and the vanishing moment 
method. The use of Newton's method for proving existence of a solution of (1.1) 
appeared in [10] combined with a method of continuity argument and more recently 
for approximation with a direct approach by finite differences in [45] for a Monge- 
Ampere equation on the torus. Wc will extend the proof of convergence of Newton's 
method of [45] in Holder spaces on bounded smooth domains. We then characterize 
the Newton's iterates as solutions of variational problems and a solution of (1.1) is 
also shown to be characterized by a variational formulation for which we derive error 
estimates for finite element approximations. 

3.1. Newton's method. We denote by C'^(O) the set of all functions having all 
derivatives of order < k continuous on Q where A; is a nonnegative integer or infinity 
and by C*^(0), the set of all functions in C^(ri) whose derivatives of order < k have 
continuous extensions to fl. The norm in C'^(O) is given by 



\u 



j=0 



A function u is said to be uniformly Holder continuous with exponent a, < a < 1 
in Q, if the quantity 

\u{x)-u{y)\ 
\x — y\'^ 

is finite. The space C'''°'{fl) consists of functions whose k-th order derivatives are 
uniformly Holder continuous with exponent a in Q. It is a Banach space with norm 

\D^u(x) — D^u{y)\ 

MW.c^m = II^^IICKS^) +SUP|/3|=feSUp^^^ \^_„\a ■ 

Next note that for any n x n matrix A, detA — l/n(cof A) : A, where coiA is the 
matrix of cof actors of A and for two n x n matrices M,N, M : N — Yllj=i ^ij^ij 
is the Kronccker product of M and N . For any sufficiently smooth matrix field 
A and vector field f, divA-^f = (divA) ■ v + A : Dv. Here the divergence of a 
matrix field is the divergence operator applied row-wise. If we put v = Du, then 
det D'^u = l/n(cof D'^u) : D^u = 1/n {coi Dv) : Dv and div(cof Dv^v — div(cof Dv)- 
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V + (cof Dv) : Dv. But divcof Dv = 0, c.f. for example [31] p. 440. Hence since D^u 
and cof D^u are symmetric matrices 

(3.4) det D'^u = -(cof D'^u) : D'^u ^ - div ((cof D^u)Du). 

n n 

Put F{u) = det D^u - f. The operator F maps C""'" into C™-2'",m > 2. This can 
be seen from the properties of the Holder norm of a product [22], p. 18. We have 
F'{u)w — (cof D'^u) : D'^w — div ((cof D'^u)Dw) for u,w sufficiently smooth. The 
Newton's iterates can then equivalently be written 

F'{uk){uk+i -Uk) 

div ((cof D'^Uk){Duk+i - Duk) 

(cof D'^Uk) : {D'^Uk+i - D'^Uk) 

We will use the last expression as it requires less regularity on Uk- More precisely, 
we will consider the following damped version of Newton's method: Given an initial 
guess uq, we consider the sequence Uk defined by 

(3.5) (cof D^Uk) : D^Ok = -(/ - fk), fk = det D'^Uk, 9k = Uk+i - Uk, 

T 

for a parameter r > 1. Our numerical results however use only r = 1. 

We recall that the domain Q is uniformly convex [40], if there exists a number rriQ > 
such that through every point Xq G dfl, there passes a supporting hyperplane tt of 
fl satisfying dist (a;,7r) > mo\x — xo\'^ ior x E dO,. We will need the following global 
regularity result, [56]. 

Theorem 3.1. Let Q be a uniformly convex domain in BJ\ with boundary in . 
Suppose (f) e C^(fi), inf f > 0, and f e C"' for some a e (0,1). Then (1.1) has a 
convex solution u which satisfies the a priori estimate 

Il'"llc2."(n) ^ 

where C depends only onn,a, inf f,fl, \ \f\\cc(U) ll^llc^- 

According to [56], all assumptions in the above theorem are sharp. We have the 
following analogue of Theorem 2.1 in [45]. 

Theorem 3.2. Let be a uniformly convex domain in B? , with boundary in . Let 
< m < f < M,f e C° for some m,M > and a G (0, 1). There exists r > 1 
depending on m,f, such that if Uk is the sequence defined by (3.5), it converges in 
C^'^ to a solution u of (1.1), for every (5 < a. 



= -Fiuk), 

= --div ((cof D\k)Duk)+f, 

Th 

= - det D^Uk + f. 



Proof. The proof essentially depends on global Holder regularity of the solution of the 
Monge- Ampere equation. Theorem (3.1) essentially gives the conditions under which 
such a regularity result holds on a bounded domain. Part of the proof has been more 
or less repeated in [35]. We note that the iterates may converge to either a concave 
or a convex solution if both exists. □ 
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3.2. Variational formulation. Using the divergence form of (3.5), the iterates can 
be characterized as solutions of the problems: Find Uk+i G H^{fl),Uk+i = g on dfl 
and 



(3.6) 



/ (cof D^Uk)Duk+i ■ Dw dx = — / (cof D^u^jDuk ■ Dw dx 

Jn n Jq 

+ [ fw dx, Vw e H^{n). 
Jn 



With an initial guess uq which solves Au = nf^^'^, for / in C", < a < 1, we have 
a sequence of uniformly elliptic problems, (see proof of Theorem 2.1 in [45]) and the 
problems (3.5) and (3.6) are equivalent. We then know that the iterates converge to 
a solution of (1.1) which solves the formal variational hmit of (3.6): Find u e i?"(Q), 
u — g on dQ such that 

(3.7) / (cof D\)Du -Dw dx = -n [ fw dx, Ww G i/"(l]) n H^in). 

The problem (3.7) is not well posed in general. For example if = 0, then both u 
and —u are solutions. 

To see that the left hand side is bounded for u e H^{fl), notice that for n = 2 



/ (cof D^u)Du ■ Dw 
Jn 



< C\\D'^u\\L2(n)\\Du\\L4(^n)\\Dw\\Li(n)- 



Next for u G H'^{Q),du/dXi G H^{Q),i = 1, . . . ,n and by the compactness of the 
embedding of H^{Q) in L'^{Q) for q > 1 when n = 2, the right hand side above is 
bounded by C||D^M||i2(Q)||M||j|^2(j^)||w||j:^2(Q). However in three dimensions, the term 
cof D'^u involves the product of two second order derivatives. For it to be bounded, 
we will need u G H^{Q) so that d'^u/dxidxj,i,j, = 1, . . . ,n G H^{Q) and we can use 
the embedding of H^{Q) in L'^{Q) for 1 < g < 6 when n = 3. We give below error 
estimates only for the two dimensional case using techniques borrowed from [32]. We 
note that in the spline element method, it is possible to impose the continuity re- 
quirements for a conforming finite element subspace of H^[Q). However for a smooth 
solution, continuity was enough for numerical evidence of convergence. 

3.3. Error estimate for 2D conforming approximation of a smooth solution. 

In this section, wc will assume that is a two dimensional domain. Put V = H'^{Q) 
and Vo = H^{Q) fl Hq{Q). Let be a conforming finite element subspace of H'^{Q), 
Vq be a conforming finite element subspace of H'^[Q) n HqIQ) with approximation 
properties 

(3.8) ini^^^vh\\v - VhWj < Cih"'^\v\\4, j = 0,1,2 

for all V G H'^IQ) for some p > 4. 

For example, in this paper, we take Vh as the spline space of degree d and smoothness 
r > 1 

^^T) = {pe C'Xn), Pit G P,, Vt G T}, 

where Pd denotes the space of polynomials of degree d in two variables and T denotes 
the triangulation of the domain Q. It is known that [43], for d > 3r+2 and < m < d. 
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there exists a linear quasi-interpolation operator Q mapping Li{Q) into the spHne 
space S^{T) and a constant C such that if / is in the Sobolev space H"^~^^{fl), 

(3.9) \f-Qf\k<Ch"'^'-'\f\m+u 

for < k < m. For our purposes, m = 3 and p = 4. If f2 is convex, the constant C 
depends only on d, m and on the smallest angle 9a in T. In the non convex case, C 
depends only on the Lipschitz constant associated with the boundary of Q. It is also 
known c.f. [30] that the full approximation property for spline spaces holds for certain 
combinations of d and r on special triangulations. We have the following theorem 

Theorem 3.3. Assume that problem (1.1) has a solution u e iJ^(Q) hence in C^(Q) 
by Sobolev embedding, then the discrete problem, find Uh & Vh, Uh — g on such 
that 

(3.10) -I [ {cof DW)Duh ■ Dwh dx^ I fwn dx, ^w^ G V^^ 

^ Jn Jo. 

has a unique solution in a neighborhood of IhU where is an interpolation operator 
associated with Vh- Moreover \ \uh — Ihiu)\\2 is at least 0(h). 

The proof of the above theorem follows the combined fixed point and linearization 
method used in [32] . The difference here is the assumption that the solution is smooth 
and the use of an inverse inequality. We start with some preliminaries. 

We consider the linear problem: Find v e Hl{Q) such that 

(3.11) / {coi D'^u)Dv Dw = [ (j)w,yw e H^{n), 
Jn Jq 

for (f) e L2(Q). 

Since u e C'^{9.),3m, M > such that 

du 

m < -^-^{x) < M, i,j = l,...,2,Va; e n. 

OXiOXj 

The existence and uniqueness of a solution to (3.11) follows immediately from Lax- 
Milgram lemma. Similarly, there exists a unique solution to the problem: Find 
Vh e Vq such that 

/ (cof D\)Dvh -Dwh^ I (t>Wh,ywh e V^. 

(,3.12j Jo, Jq 

f = on dQ. 

We note that the constant C above depends on u and that since Q is assumed convex, 
V G H'^{Q) by elliptic regularity. 

We define a bilinear form on -f^o(^^) x Hq{Q) by 

(3.13) B[v, w]^ (cof D\)Dv ■ Dw, 

Jn 

and for any Vh G V^, Vh = g on dQ, we define T{vh) by 

(3.14) B[vh - T{vh),Wh] [ (cof D%h)Dvh ■ D^^ dx + f f^h dx, VV'^ e V^K 
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Since the solution of the above problems exists and is in V^. T{vh) G V^^ T{vh) = g 
on dfl. A fixed point of the nonlinear operator T corresponds to a solution of (3.10) 
and conversely if Vh is a solution of (3.10), then Vh is a fixed point of T. We will show 
that T has a unique fixed point in a neighborhood of Ih{u). Put 

Bh{p) = {vh e Vh.Vh ^goTidVL,\\vh - 4ii||2 < p}- 
First, we note that 

(3.15) B[vn-T{vh),Wh] = - [ det V^;, dx + [ f^n dx, y^/j^ G 
Lemma 3.4. There exists C2 > such that 

\\hu-T{hu)\\2<C2hP~^\\u\U. 

Proof. Put Wh — IhU — T{Ihu). We have using det D^u — /, 

B[wfi,Wfi]— / (^det D^u — det I hujwfidx. 
Jn 

Now, let Il{u) be a moUifier of IhU. We have 

B[wh, Wh]^ { det D\ - det D'^Il{u))wh dx + { det D^Il{u) - det D'^Ihu)wh dx 
Jq Jn 

= j (^{coitD'^u+{l-t)D^Il{u)) : D'^{u- Il{u))^Whdx 

+ / (det D'^Il{u) - det D^Ihu)wh dx for some t e [0, 1] 

< I |cof tD\ + (1 - t)D^Il{u))\U\u - Il{u) MwhWo 

+ 1 1 det D''ll{u) - det D'^huMwhWo 

< C\\D'^u\\oo\\u - h{u)\\2\\wh\\o < CM/i^'-^||M||4||w^||o, 

since || det D'^If^{u) — deiD'^Ihu\\Q — >■ as e — >■ 0. We conclude that 

C 

\\wh\\\ < ChP'^\\u\\4\\wh\\o, and \\wh\\2 < -^\\wh\\i < Ch^~^\\u\\4, 

using the cocrcivity of the bilinear form B with a constant C which depends on m 
and M and an inverse estimate. □ 



Lemma 3.5. There exists Hq > andO < po(^o) such thatT is a contraction mapping 
in the ball Bh{po) with a contraction factor 1/2. 
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Proof. For Vh,Wh £ Bh{po), and iph G ^o, let vl and w\ denotes moUifiers for Vh and 
Wfi respectively. 

B[T{vh) -T{wh),iph] ^ B[T{vh) - Vh,iph] + B[vh - Wh,iJh] + B[wh - T{wh),iph] 

(det D^Vh — det D'^Whji^h dx 

n 

+ 



/ (cof D\){Dvh - Dwh) ■ Di/jh dx 

Ae+ {coi D^u){Dvh - Dwh) ■ Diph dx 
Jn 

+ / (det D^v'^ - det D^wDiph dx, 



where = J^{det D'^Vh — det D'^vl)ijjhdx + J^{det D'^Wh — det D'^wl)ijjh dx ^ as 
e ^ 0. We have for some t e [0, 1]), 

B[T{vh) - T{wh),i^h\ =A,+ [ (cof D\){Dvh - Dwn) ■ D^Jh dx 

Jn 

- [ (cof tDX + (1 - t)D'wI){DvI - DwDi^n dx 
Jn 

= / (cof D\ - (cof tD'^vl + (1 - T)D^wl)){Dvh - Dwh) ■ Diph dx 
Jn 

+ A, + B„ 

where B, = ^(cof tD%l + (1 - T)D^wl){Dvh - Dwh - Dvl + DwDipndx ^ as 
e ^ 0. Put = cof D'^u - (cof tD'^v^ + (1 - t)D'^wI). We have 

1 1 lo = \\D^u - D'vl + T{D'wl - DX) I lo < I \D'u - D^l lo + I \D'wl - D^W, 

<\\D^u- D^huWo + WD'lhU - D\,\\o + I \D\h - Z^Xl lo + I I^H - D^^JhWo 

+ \\D^WH - DWWo + WDh'h D^jiWo 
< ChP~^\\u\\4 + 3po + 2\\vh -vl\\2 + \\wh- <| I2. 

As e ^ 0, we obtain 

B[T{vh)-T{wh),M<C{hP-^\\u\U + Po)\\vh-WhM\i^h\m- 
By coercivity and an inverse estimate, 

I \T{v^) - T{wh) 1 12 < C{N>-^\\u\ U + ^)\\VH-Wh\\2. 

First choose h such that /i^~^||M||4 < 1/4 then po ^ h/A. The result follows □ 

Proof, (of Theorem 3.3) Let pi — 2C2hP''^\\u\\i. We first show that T maps Bh{pi) 
into itself. For Vh £ Bfi{pi), 

\\hu-T{vh)\\2<\\IhU-T{Inu)\\2 + \\T{hu)-T{vh)\\2 < ^ + l\\hu - vj^lU 
^ Pi .Pi 



11 

By the Brouwer's fixed point theorem, T has a unique fixed point in Bh[pi) which is 
Uh, i.e. T{uh) = Uh. Next, 

\\U - Uh\\2 < 11^^ - hu\\2 + Whu - Uh\\2 < CihP~'^\\u\\4 + \ \IhU - T{Uh)\\2 

< CihP-^\\u\\4 + 2C2hP-^\\u\\4 < ChP-^\\u\\^, 
for h sufficiently small. We have the result. □ 

3.4. Vanishing moment methodology. The vanishing moment methodology ap- 
proach to (1.1), consists in computing a solution of the singular perturbation problem 

(3.16) -eA\ + det D^u = /, in Q, u = g, Au = on dQ. 

It is an analogue of a singular perturbation problem 

oil 

eA\ - Aii = / in Q, = 0, — = 0, on dn 

on 

which was addressed numerically in [47] and also in the spline element method [6]. 
The analogy holds as many properties of the Laplace operator have a counterpart for 
the Monge- Ampere operator. 

The Newton's iterates in the vanishing moment formulation consisting in solving the 
sequence of problems: Find u^+i € H"-{fl) with u^+i = g on dfl 

(3.17) 

/ Auk+iAv dx + / (cof D^Uk)Duk+\ ■ Dv dx — / (cof D^Uk)Duk ■ Dv dx 

Jn Jn ^ Jn 



e 



+6^ f ^ds- f fv dx, Vv e H^{n) n H^{n). 

J an Jn 



Formally as e approaches 0, the solution of the above problem degenerates to the 
solution of (3.6), a result which will be illustrated numerically in the next section. 



4. Numerical results 

The first two subsections are devoted to two dimensional and three dimensional nu- 
merical results respectively The three methods are compared on three test functions 
for 2D experiments. 

Test 1: A smooth solution u{x,y) = e^^^^+^^^Z^ so that f{x,y) = {1 + x'^ + y^jgCaj^+y^) 
and g{x,y) = 6*^^^+^^)/^ on dQ. 

Test 2: A solution not in H'^(Q.), u{x,y) = — -^2 — x"^ — y"^ so that f{x,y) = 2/(2 — 

2;2 _ y2-^2 g^^^ y-j ^ — 1^2 — — on dQ. 

Test 3: No exact solution is known. Solutions are either convex or concave. Here 
f{x,y) = 1 and g{x,y) = 0. 

We give numerical evidence of the robustness of the spline element method for the 
singular perturbation problem associated to the vanishing moment methodology. For- 
mally as e approaches 0, the problem (3.16) degenerates to the problem (1.1), which 
can be solved by Newton's method when a smooth solution exists. We show here 
numerically that the solution of (3.17) converges to that of (3.6) as e approaches 0. 
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Unlike [34], there is no boundary layers issue with the spline element method. These 
results are illustrated in Tables 4.1 and 4.2. 

In general, vanishing moment (with Newton) gives Newton's method result for e ~ 
10~^. In particular, Newton's method diverges for the non smooth solutions of Test 
2 and 3. However with e large, which implies that the equation solved is much 
further from the actual problem, divergence can be avoided in the vanishing moment 
methodology. We refer to the problem (1.1) as reduced in the tables. 

In the two dimensional case of Test 3, both concave and convex solutions can be 
computed by either changing the initial guess or the structure of the approximations. 

(1) Newton's method: initial guess iwo) 

(2) Vanishing moment: ±e and initial guess ±-Uo, 

(3) BFO iterative method: Uk+i = ±^/{Auk)'^ + 2(/ - dct D'^Uk). 

Since Newton's method diverges for Test 3, we illustrate this feature of the method on 
a fourth test on a non-square domain. This also helps contrast with finite difference 
methods. 

Test 4: The domain is the unit circle which is discretized with a Delanauy triangu- 
lation with 824 triangles and the test functions are u{x, y) — e^^ +^ (convex) and 
u{x, y) = —e^^ "'"^ (concave) with the corresponding right hand side and boundary 
conditions. 

Since none of the methods perform convincingly on Test 2 in the spline element 
framework, the methods are tested for the three dimensional case on two other test 
functions analogues of Test 1 and Test 3. We only consider the performance of the 
BFO and vanishing moment method. 

Test 5: u{x, y, z) = e(^'+2''+^')/=^ so that /(x, y, z) = 8/81(3+2(x2+yH^^))e(^'+2''+^') 
and g{x, y, z) = e^^^+y^+^^)l^ on 9Q. 

Test 6: /(x, y,z) — 1 and g{x, y, z) — 0. 

The initial guess of the Newton's iterations is the solution of the Poisson equation 
Am = n/^/" n = 2, 3 in r2, u — gon dQ. 

We also illustrate the performance of the 3D BFO method on a degenerate Monge- 
Ampere equation. 

Test 7: f{x, y,z) = and g{x, y, z) = \x — l/2\ . 

In the two dimensional case, to approximate a concave solution, one should solve 
Au ~ —2\fj. But unless m = on 9Q, as in Test 3, it is not clear which boundary 
condition to use. Note that if m is a smooth convex function. Aw > 0. To compute 
the concave solution of Test 4, we first solved the Monge-Ampere equation with 
the negative of the boundary condition, then use the negative of that solution as an 
initial guess. However, a good initial guess could not be found if we uniformly refine a 
Delanauy triangulation of the circle with 143 triangles, but convergence was obtained 
with the choice in Test 4, perhaps because the domain is closer to being smooth. 
For the vanishing moment calculations, the initial guess was taken as the biharmonic 
regularization of a suitable Poisson equation, for example, — eA^w + Am = n/^/" n = 
2, 3 in Q, u — g, Au — on dO,. We simply took the zero function as initial guess 
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d 


L'^ norm 


norm 


norm 


3 


1.0610 10-^ 


1.1101 10-^ 


1.6383 10-^ 


4 


3.5127 10"^ 


4.8553 10-4 


9.0596 10-^ 


5 


4.1572 IQ-*^ 


6.5142 IQ-^ 


1.9364 10-^ 


6 


1.9685 10"'^ 


3.6401 10-^ 


1.4774 10-4 


7 


2.2699 10-« 


4.1498 10-'' 


2.2424 10-^ 


8 


1.2430 10-y 


2.2586 10-« 


1.5479 10-*^ 



Table 1. 2D Newton's method for Test 1, /i = 1/2 



d 


norm 


norm 


H'^ norm 


3 


1.2809 10-4 


2.6554 10-^ 


8.9587 10-^ 


4 


1.6278 10-^ 


4.5619 10-^ 


1.7395 10-^ 


5 


1.1531 10"'' 


2.3916 10-*^ 


1.3444 10-4 


6 


1.7705 10-y 


6.8506 10-« 


5.5403 10-*^ 


7 


1.4548 10-1° 


3.7545 10-9 


3.9490 10-^ 


8 


8.1014 10-1^ 


5.3353 10-1" 


7.2159 10-« 



Table 2. 2D Newton's method for Test 1, /i = 1/4 

in the BFO method. Unless otherwise stated, we use splines for all the numerical 
experiments. Even for the BFO iterative method which requires only solving Poisson 
equations as in that case we obtained better numerical results (smooth graphs) for 
Test 3. We hsted Uit, the number of iterations of the BFO method. We do not claim 
that our numerical solutions are convex but that they approximate convex functions. 
Convexity (or concavity) is not explicitly enforced in the numerical iterations. 

4.1. Two-dimensional Monge- Ampere equation. The computational domain is 
the unit square [0, 1]^ which is first divided into squares of side length h. Then each 
square is divided into two triangles by the diagonal with negative slope. As evidenced 
in the last fine of Table 4.1, we noted a degradation of the performance of the BFO 
iterative method even for a smooth solution when the number of degrees of freedom 
is large, either for smaller h or higher degree. This may be an indication that the 
method is not suitable for a general finite element implementation but is more likely 
a consequence of the conditioning properties of the iterative method (2.2). For Test 
2, we display the results for cubic splines, but much higher order accuracy, of the 
order of 10-^ was obtained with C° splines. We caution that in our implementation, 
this may lead to non smooth numerical solutions. 

For Test 3, we displayed both the graph and the contour of both convex and concave 
solutions. To get good results with the vanishing moment method, we experimented 
with a combination of the parameters e and h. 

4.2. Three-dimensional Monge- Ampere equation. We used two computational 
domains both on the unit cube [0,1]^ which is first divided into six tetrahedra (Domain 
1 for Test 4) or twelve tetrahedra (Domain 2 for Test 5) forming a tetrahedral partition 
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e 


norm 


norm 


norm 


10"^ 


1.727110--^ 


2.091010-^ 


8.633810-1 


10-4 


1.856310-4 


3.558110-^ 


2.005010-1 


10-^ 


1.891710-^ 


4.070010-4 


2.411910-^ 


lO"'* 


1.818210-^ 


4.077510"^ 


2.438810-^ 


10-'^ 


1.244110-'^ 


4.095110-^ 


2.494910-4 


io-« 


1.011910-'^ 


2.296210-^' 


1.302910-4 


10-9 


1.138410-^ 


2.379010-6 


1.338210-4 


10-10 


1.151610-^ 


2.390310-'^ 


1.343810-4 


10-11 


1.153010"^ 


2.391410-6 


1.344310-4 


10-14 


1.153110-'' 


2.391610-6 


1.344410-4 


Reduced 


1.1531 lO"*" 


2.3916 10-6 


1.3444 10-4 



Table 3. 2D numerical robustness Test 1, h — d — h. 



h 


nit 


norm 


iJi norm 


norm 


1/2^ 


41 


2.8275 10-6 


6.1372 10-5 


1.8845 10-3 


1/2^ 


37 


5.4642 10-« 


2.1971 10-6 


1.2972 10-4 


1/2^ 


38 


8.3164 10-1" 


7.2252 10-« 


8.4790 10-6 


1/2^ 


37 


2.7871 lO-'^ 


1.4089 10-« 


1.0809 10-6 



Table 4. BFO iterative method for Test 1, d = 5 



h 


norm 


norm 


l/2i 


2.195410-^ 


1.640910-1 


1/2^ 


3.609710-'^ 


6.140510-^ 


1/2^ 


1.068510-3 


4.097810-^ 


1/24 


5.083810-3 


2.804810-1 


1/2^ 


2.579710+3 


2.268810+5 


1/26 


1.845210+4 


3.592210+6 



h 




norm 


ifi norm 


l/2i 


50 


2.3921 10-1 


1.1900 


1/2^ 


159 


1.2585 10-1 


7.1292 10-1 


1/23 


151 


1.0341 10-1 


6.4299 10-1 


1/24 


160 


9.6031 10-2 


6.2088 10-1 


1/2^ 


199 


9.4551 10-2 


6.2453 10-1 


1/26 


8 


1.6977 10-2 


2.2925 10-1 



Table 5. Newton's method and BFO iterative method for Test 2, d — ?> 



h 


L2 norm 


iJi norm 


l/2i 


7.668010-3 


7.449110-2 


1/22 


1.453610-3 


3.924410-2 


1/23 


9.872710-3 


2.511210-1 


l/2'i 


5.681910-3 


2.492710-1 


1/2^ 


1.9830 10+4 


1.1812 10+6 



h 


L2 norm 


iJi norm 


l/2i 


7.825410-3 


9.318410-2 


1/22 


1.064610-2 


9.520110-2 


1/23 


1.130610-2 


9.615410-2 


1/2^ 


1.150010-2 


9.133610-2 


1/25 


1.162510-2 


8.778510-2 


1/26 


1.168110-2 


8.563210-2 



Table 6. Vanishing moment Test 2 e = 10-^ and e = 10-2, d = 5 
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Figure 1. Vanishing moment on Test 3, h = 1/2^, d = 3,e = 10~^. 




Figure 2. Vanishing moment on Test 3, h = 1/2^, = 5,e = -10"^. 




Figure 3. BFO iterative method on Test 3, h = 1/2"^, d = 3. 



Ti- This partition is uniformly refined following a strategy introduced in [4] similar 
to the one of [51] resulting in successive level of refinements Tk, k = 2,3, . . .. For Test 
6, we plot the graph of the function as well as its contour in the plane x = 1/2 as 
well as slices in the x— direction. 
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Figure 4. BFO iterative method on Test 3, h = 1/2"^, d = 3. 




Figure 5. Approximations of smooth concave and convex solutions 
on a non rectangular domain, Test 4, = 5, r = 1 



d 


norm 


norm 


norm 


3 


1.2338 10-^ 


7.6984 10-^ 


4.4411 10-^ 


4 


1.6289 10-3 


1.4719 10-2 


1.3983 10-1 


5 


1.5333 10-3 


8.7312 10-3 


6.0412 10-2 


6 


1.2324 10-^ 


9.7171 10-* 


1.0584 10-2 



Table 7. Newton's method Test 5, Domain 1 on Xi 



d 


norm 


norm 


norm 


3 


3.1739 10-3 


2.3005 10-2 


2.4496 10-1 


4 


3.2786 10-* 


3.5626 10-3 


5.2079 10-2 


5 


2.4027 10-^ 


3.9210 10-* 


8.8868 10-3 


6 


1.3821 10-*^ 


2.2369 10-^ 


6.0918 10-* 



Table 8. Newton's method Test 5, Domain 1 on 72 




Figure 8. Slices in the x— direction Test 6 on Domain 2 and X3, d = 3 
, Vanishing moment r = 2, e = 10~^ and BFO d = 5,r = 1 
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e 


Li^ norm 


norm 


H'^ norm 


10-^ 


6.6870 10"^ 


3.9292 10-1 


2.8852 


10-^ 


1.8832 10"^ 


1.3137 10-1 


1.5882 


10^^ 


2.4237 10^^ 


2.5273 10-^ 


5.3206 10-1 


10^4 


2.5661 10-4 


3.2633 10-^ 


7.9936 10-^ 


10-^ 


3.1058 10-^ 


5.0367 10-4 


1.2543 10-^ 


10-^ 


2.3519 10-^ 


3.9165 10-4 


8.9744 10-^ 


10-^ 


2.3964 10-5 


3.9193 10-4 


8.8921 10-3 




2.4027 10-5 


3.9210 10-4 


8.8868 10-3 


Reduced 


2.4027 10-5 


3.9210 10-4 


8.8868 10-3 



Table 9. 3D numerical robustness Test 5, Domain 1 on 72, li = 5 



2 


3.1739 10-3 


2.3005 10-^ 


2.4496 10-1 


3 


1.6859 10-^ 


1.0519 10-1 


9.1615 10-1 


59 


1.1283 10-3 


7.1385 10-3 


7.3671 10-^ 


38 


2.1423 10-4 


1.4452 10-3 


1.8083 10-^ 


35 


4.5582 10-5 


3.0440 10-4 


4.0506 10-3 



Table 10. BFO iterative method for Test 5, on 72, ci = 5 




Figure 9. BFO Test 7 on Domain 2 and X3, d = 5,r = 1 

5. Concluding Remarks 

Remark 5.1. For the finite element approximation of (1.1), we note the remark in 
[27], "Newton's and conjugate gradients methods may he well- suited for the solution 
of . . . combines the difficulty of both harmonic and bi-harmonic problems, making the 
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approximation a delicate matter, albeit solvable . . . If . . . has no solution, we can ex- 
pect the divergence of the Newton ..." We have established that Newton's method 
performed well for smooth solutions. Another problem which also combines the diffi- 
culty of both the harmonic and biharmonic problem is a singular perturbation problem 
we addressed in [6] by the spline element method. Here it is seen that for the singular 
perturbation problem arising from the vanishing moment methodology, the spline ele- 
ment method is robust. Moreover we note that numerical results for smooth solutions 
using Newton's method in the spline element method are more accurate than what can 
be achieved using Argyris elements and the vanishing moment methodology [34]. 

Remcirk 5.2. It is still not known whether the BFO iterative method always con- 
verges even in the case of smooth solution. Nor is known whether Newton's method 
always converges on a non-smooth domain. We have not addressed the convergence 
of the vanishing moment methodology to viscosity solutions as these results have been 
announced in [34] . 

Remark 5.3. We used cubic splines on most of the approximations even though 
they do not have full approximation power on general meshes. This reduced the com- 
putational cost. One may use for full approximation power, special meshes as in [30] 
or [53]. 

Remark 5.4. The BFO iterative method introduced in [16] was tested on some very 
singular right hand sides. It was noted that it is slower than another iterative method 
based on a different algebraic manipulation of the Monge- Ampere equation. The latter 
does not seem directly amenable to finite element computations. We note that for 
singular sources, specialized finite elements may have to be used and even in the finite 
difference context, specialized finite difference methods [55] or fast Poisson solvers [54] 
or precondition ers could have been used. 

Remcirk 5.5. The problem: Find u such that 

(5.1) I {cofD'^u)Du-Dv dx^-n I fv dx, 

Jn Jq 

is the Euler- Lagrange equation of the functional 

(5.2) J{v)^ I icofD'^v)DvDvdx + 2n / fvdx. 

Jq Jq 

If V — on dfl, we have 

J[v) ^-n I (det D'^v)v + 2n fvdx, 
Jq Jq 

and a generalized solution of (1.1) has been shown in [11, 57] to be a minimizer of a 

related functional on the set of convex functions. 
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